The idea is very simple. Start with a unit vector pointing at the north pole (1,0,...,0). Then generate random numbers from a standard normal and scale them so that they have a unit length. To put it differently, a sample of n values from the uniform distribution on the sphere is generated. Then calculate the rotation matrix required to go from the north pole to each of a generated vector.
References
Amaral G.J.A., Dryden I.L. and Wood A.T.A. (2007).
Pivotal Bootstrap Methods for k-Sample Problems in Directional Statistics and Shape Analysis.
Journal of the American Statistical Association, 102(478): 695--707.